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We analyse the buckling stability of a thin, viscous sheet when subject to simple shear, 
providing conditions for the onset of the dominant out-of-plane modes using two mod- 
els: (i) an asymptotic theory for the dynamics of a viscous plate and (ii) the full Stokes 
equations. In either case, the plate is stabilised by a combination of viscous resistance, 
surface tension and buoyancy relative to an underlying denser fluid. In the limit of vanish- 
ing thickness, plates buckle at a shear rate 7/ (/id) independent of buoyancy, where 2d is 
the plate thickness, 7 is the average surface tension between the upper and lower surfaces 
and \i is the fluid viscosity. For thicker plates stabilised by an equal surface tension at 
the upper and lower surfaces, at and above onset, the most unstable mode has moderate 
wavelength, is stationary in the frame of the centre-line, spans the width of the plate with 
crests and troughs aligned at approximately 45° to the walls and closely resembles elastic 
shear modes. The thickest plates that can buckle have an aspect ratio (thickness/ width) 
approximately 0.6 and are stabilised only by internal viscous resistance. We show that 
the viscous plate model can only accurately describe the onset of buckling for vanish- 
ingly thin plates but provides an excellent description of the most unstable mode above 
onset. Finally, we show that by modifying the plate model to incorporate advection and 
make the model material frame- invariant, it is possible to extend its predictive power to 
describe relatively short, travelling waves. 



1. Introduction 

Folding, buckling and coiling are phenomena frequently associated with thin elastic 
solids. However they also occur in very viscous films and filaments whenever compression 
is faster than can be accommodated by film or filament thickening. Viscous buckling has 
been studied in a variety of contexts over the last half century. A primary motivation 
for some of the earliest work was understanding the buckling of layered geological strata 
modeled as very viscous fluid layers (w ith viscosit ies that range from 10 16 to 10 21 Pas) 
This work was pioneered by Biot (e.g.. [Biotlll96lh . who examined the two-dimensional, 
small-deformation folding of viscous layers embedded in a less viscous medium and sub- 
jected to layer-parallel compression. He used t he Stokes-Ra yleigh analogy relating viscous 
creeping flows with their clastic counterparts (IStruttlll945h and the concomitant similar- 



ity between elastic and viscous governing equations to develop expressio ns for the criti 



cal load and w avelength of the instability. Many subsequent studies (e.q, lR ambcr gll963l; 
Chapplel 1968 ) added further physical effects (a summary is given bv I Johnson fc Fletcher 



1994h . Viscous buckling is also encountered in more familiar contexts: the folding of cake 
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batter pouring into a pan, the wrinkling of a layer of cream on hot milk, and the coiling 
of a stream of honey falling from a spoon. Each of these examples also have industrial 
analogues in the spinning of polymeric fibres and in the shaping and blowing of glass 
sheets and shells. This second set of applications has provided a new impetus to the 
study of these problems using a combination of approaches. 

At a theoretical level, a systematic asymptotic reduct ion of the full governing equa- 
tions to the thin geometry of interest was carried out by iBuckmaster. Nachman fe Ting 
(|1975l ). who investigat ed the l arge amplitu de deformation of a filament. Building on 
their scaling relations, iHowelll (|1994 Il996l ) developed asymptotic equations governing 
the evolution of thin filaments and sheets in a variety of scenarios. In particular he de- 
rived equations for small deformations of viscous sheets equivalent to the Foppl-von 
Karman eq uations for elastic plates (a linearized version thereof had been stated by 
analogy by IBeniamin fe Mullinl [1988|). Subsequently, a number of bending, stretching 
and buckling phenomena inv olving viscous filaments and planar deformations of viscous 
sheets have been explained: Yarin fe Tchavdarovl (1996) conside r ed the onset of buck- 
ling in a filament impinging on a wall, iTeichman fe Mahadevan (2003) considered the 
viscous catenary using a combinat ion of scaling, asymptotic and numerical approaches; 
Mahadevan. Ryu fe Samuel (1998) and ISkorobogativ fc Maha devanl (120001) provided a 



simple physical picture for the the dif ferent regimes of coiling and fo lding of filaments on 
impact with a stationary surface; and IChiu- Webster fc Listen (|2006l ) considered the com- 
plex 'stitching' patterns of a filament impacting a moving s urface. Somewhat fewer studies 
have i nvestigated three-dimensional deformation of sheets: Silveira. Chaieb fe Mahadevanl 
(|2000h considered t he wr inkling of a ruptured viscous bubble collapsing under its own 
weight; Teichman] ( 2002 ) consider ed th e buck ling of sheared viscous sheets in both a 
rectilinear and Couette geometry; Ribe (|200ll ) derived asymptotic equations for sheets 
of high curvature and analysed aspects of geophysical problems such as trench roll back; 
Slim et all (l2009h briefly co nsidered buckling of a thin vis c ous s heet by an underlying, 
less viscous fluid flow; and iMahadevan. Bendick fc Liang ( 2011 ) analysed the form of 
tectonic subduction zones. This summary is by no means comprehensive but highlights 
the evolution of, and recent interest in, viscous buckling problems, especially involving 
two dimensional deformation of a viscous plate or shell. 

Here we study the shearing of a thin, very viscous sheet in a plane Couette geometry. 
Specifically, we consider an initially uniform, thin layer of viscous Newtonian fluid of 
finite width and infinite length sheared by the constant-velocity motion of bounding 
walls. The layer floats on a deep lower fluid, which contributes intcrfacial tension and 
a gravitational restoring force. The upper surface is open to the atmosphere and only 
experiences surface tension. Contrary to the situation for an infinitely thick sheet, which 
is linearly stable to shear for all values of the shear rate, the thin sheet can and docs 
respond by buckling when sheared. We present the conditions for the onset of this linear 
in stability, as well a s growth rates and mode profiles above onset, expanding on the work 
of lTeichmar] (|2002[) . 

There is a superficial similarity between the plane Couette problem treated here and 
the circular Couette probl em of the ann ular shearing of a thin viscous film, a problem first 
studied experimentally by Tavlor (119691) . and sub sequently by others (jSuleiman fc Munson 
19811: IBeniamin fc Mullinl Il98a 



Teichman! 120021 ) . However, a fundamental difference is 



that the annular geometry naturally introduces two length scales, one associated with 
the gap and the other with curvature, while the rectilinear problem has just a single 
length scale. This difference is manifest in the pla te model predic ting a self-consistent 
onset at moderate wavelength in the annular case ( Teichmanll2002 ) and an inconsistent 
onset at infinitesimal wavelength in the rectangular case. 
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We use both the full Stokes equations and the viscous plate model to investigate growth 
rates of infinitesimal perturbations to the simply-sheared planar base state. Using the 
former, we show numerically that the thickest plate that can buckle has aspect ratio 
approximately 0.6 and is stabilised by internal viscous resistance alone. For thinner sheets 
with surface tension but no buoyancy, we establish the dependence of the critical wall 
speed on the plate width, thickness and surface tension coefficients. The viscous plate 
model is unable to reproduce these onset conditions except in the limit of vanishing plate 
thickness. Nevertheless, for plates of aspect ratio up to around 0.04, given the system 
parameters at onset, the most unstable mode for this model accurately reproduces the 
mode profile and wavelength predicted using the full Stokes equations. Above onset, 
the most unstable mode has moderate wavelength, is stationary in the frame of the 
centre-line, has crests and troughs aligned at approximately 45° to the sidewalks, and 
close l y resembles the modes of the e lastic analogue problem (e.g., Southwell k, Skaiil 



1924t iBalmforth. Craster fc S lim 2008). Waves shorter than order the plate thickness are 



suppressed by internal viscous resistance. The shortest unstable modes form a pair of 
travelling waves, each concentrated in one half of the plate and propagating at a fraction 
of the corresponding wall speed. These modes cannot be captured by the asymptotic 
viscous plate model which do not contain the advection term responsible for symmetry 
breaking. By modifying the model to include this term, we show that we can accurately 
recover the critical wavenumber and associated growth rates, propagation speeds and 
mode profiles. 

The structure of the paper is as follows. In ^5] we describe the geometry and important 
non-dimensional parameters. In $3]wc formulate the low-dimensional viscous plate model. 
Using this model in 21 we briefly describe details of pure compression to set the scene for 
the shear instability, and present a parameter space investigation of the onset of shear- 
induced buckling, as well as the growth rates and mode profiles. In f|5l we turn to the full 
Stokes description. We present the linearised perturbation equations about the simple 
shear base state, present a numerical investigation of the parameter space and compare 
our results to those of the plate model, the plate model incorpora ting advection and a 
short-wavelength approximation due to I Benjamin &: Mullinl (|l988l ). In g6] we summarise 
the parameter space before presenting our conclusions in Sj7] 



2. Geometry 

We start with a description of the geometry and the fundamental non-dimensional 
parameters. The configuration is sketched in figure [1] a thin layer ("plate") of very 
viscous, Newtonian, incompressible fluid of viscosity p and density p floats on a deep 
layer of fluid with density p~ > p. We assume inertia in the plate and viscosity in 
the underlying fluid are both negligible (Reynolds numbers in the plate are order 10 -3 
and viscosity ratios between the underlying fluid and the plate are order 10 -5 in typical 
experiments). The upper surface is open to the atmosphere. The viscous plate has initially 
uniform thickness 2d, width 2L and infinite length. We use Cartesian coordinates to 
describe the system, with the origin on the undeformed centre-line, x-axis directed along 
its length, y-axis its width and z-axis perpendicular to its undeformed centre-plane. Along 
the lateral edges y = ±L, the plate is clamped to and sheared by the bounding walls 
which move parallel to their length with velocity ±[7. 

Three non-dimensional parameters arise in the problem naturally: the aspect ratio of 
the plate, 

a = d/L, 




the scaled, inverse capillary numbers for the upper and lower surfaces and their mean: 

r± = 1 ± /(a l iu) 1 r = (r+ + r-)/2; 

where 7 ± are the coefficients of surface tension at the two surfaces, and the gravity 
numbers 

G = P gL 2 /(anU), G~~ = p~ gL 2 /(a/it/), 

where g is gravity. These measure the importance of gravity on the plate and on the 
underlying fluid respectively, relative to viscous shear in the plate. The appearance of 
the aspect ratio in the inverse capillary and gravity numbers ensures that the stabilising 
effects of surface tension and gravity scale in the same way as the destabilising effect of 
shear with variations in the plate thickness. 



3. Low-dimensional viscous plate theory: formulation 

The viscous plate equations are valid for small deflections of a very viscous fluid sheet 
whose thickness is much smaller than any extrinsic horizontal length-scale such as the 
channel width or intrinsic length-scale such as the wrinkle wavelength. For such sheets, 
out-o f-plane sinuou s deformations occur much quicker than varicose thickening and thin- 
ning (lHowelll[i~996|) . and only the former are captured by the model. Thus the sheet 



thickness remains constant at 2d, and the dynamics are best described in terms of the 
mid-plane displacement H from z = 0. A physically motivated asymptotic derivation of 
the governing equations is provided in Appendix [A] here we provide a summary. 
Balance of forces in the plane of the sheet leads to (see (|A 4[) in Appendix [AJ 

VZ = 0, (3.1) 



where the gradient operator involves only the in-plane components (x,y) (this shall be 
our convention throughout, unless explicitly stated otherwise). Here T. is the tensor of 
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t = Apd E + tr(E)l , (3.2) 

where I is the two-dimensional identity and E is the in-plane deformation rate tensor 
given by (see (|A3|) in Appendix |A"| . 

E = - (Vu h + Vfi£) + i (VHVw + VwVtf) , (3.3) 

where superscript T denotes transpose, \ih is the in-plane velocity field in the mid-surface 
of the plate and w is the out-of-plane velocity, given by 

dH , , 

w = -jr. (3.4) 
at 



We note that the pressure appears indirectly in the expression (|3.2p via a Trouton ratio 
(of four in two dimensions) and the trace of the deformation rate tensor. In (|3.3|) . the first 
two terms are due to in-plane velocity gradients while the last two describe the stretching 
rate of the mid-surface due to out-of-plane deformation and arise from differentiating 
the term (VH) 2 . Equation (|3.4j) couples the centre-plane deflection directly to the fluid 
velocity perpendicular to the sheet (equation (|A ip in Appendix [A")) . 
The vertical force balance equation yields (see (|A 5|) in Appendix [AJ 

^d 3 V 4 u> = V ■ (Z ■ VH) + ( 7 + + 7 -)V 2 H - p-gH, (3.5) 

Here the left-hand side is the Laplacian of the rate of change of mean curvature, describing 
the time-dependent resistance to bending. In conjunction with (|3.4p , it can be shown that 
this term regulates growth or decay rates of out-of-plane modes; it cannot control whether 
the system is stable or unstable. Its effect is largest for short waves. The first term on 
the right-hand side is an anisotropic Laplace pressure encapsulating the projection of the 
in-plane stresses in the out-of-plane direction; it is stabilising if the principal in-plane 
stresses are tensile and destabilising if they are compressive. The final two terms on the 
right are the stabilising effects of surface tension and buoyancy respectively. Both are 
active at all length-scales, but the former is most prominent at short wavelengths while 
the latter is most significant at long wavelengths. 

We note that the model ignores contributions due to advection, an omission which 
has two significant implications. First, the model is not material-frame invariant relative 
to translation and rotations in the plane. This is asymptotically correct in the limit 
a — > 0, provided we use a frame of reference in which the advection of perturbations 
into a region is insignificant compared to the generation of perturbations by the out-of- 
plane velocity. Second, we shall see in <j5] that the advective terms are fundamental for 
describing certain qualitative features at moderately short wavelengths. Reincorporating 
the advective terms at leading order eliminates the apparent inconsistency and extends 
the predictive power of the model. We shall describe this modification in S}5] 

Along the lateral boundaries 

u/i = (±C/,0), w = dw/dy = 7 on y = ±L. (3.6a,b,c) 
3.1. Scaled equations 

We scale using the sheet half-width I as a length and the wall speed U as a velocity, and 
include factors of the aspect ratio appropriate for a thin sheet (see Appendix [A} using 
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hats to denote dimensional variables. Thus we set 

u h = Uu h , w = (U/a)w, (x, y) = L(x, y), H = aLH, 

E = (U/L)E, t = anlTL, t = a 2 {L/U)t, (3.7) 
reducing the system (|3.ip - (|3.6[) to 



o r) fj 

VE = 0, -V 4 w = V-(E- VH) + 2TV 2 H -G-H, — = w, (3.8a,b,c) 

O Ob 

E = 4[E + tr(E)l], E = ^ (Vu h + Vu[ + VHVw + VwVH) , (3.8d,e) 
with boundary conditions 

u h = (±1, 0), w = dw/dy = on y = ±1. (3.8/,.g 7 /i) 
The non-dimensional surface tension and buoyancy parameters T and G - are as defined 

in m 

3.2. Analogy with an elastic plate 

There is a close connection between the gover ning equation s for elastic and viscous plates, 
following from the Stokes- Ray leigh analogy (|Struttlll945l ). and the associated buckling 
instabilities that they describe, which we now discuss. For an incompressible elastic ma- 
terial wath_Yo4mg^sjnodjrhrs_F^^ Karman equations are given by 
(e.g., Timoshenko k, Woinowskv-Krieger 19591 ) 

V • E = 0, ^Yd 3 \7 4 H = V • (E • VH) + T S7 2 H - p~gH, 



E + tr(E)l 



1 



Vu h + Vui + VHVH 



where we now associate E with the in-plane deformation tensor and ii/j with the in- 
plane displacement of the centre-plane. To complete the analogy with the viscous plate 
described above, we have included T , an isotropic background tension (T > 0) or 
compression (To < 0), and a buoyant restoring force. This system becomes identical to 
the viscous plate model on identifying \idjdt with Y/3 and To with 7 + +7~. Boundary 
conditions (|3.6p translate into clamped edges 



Q h = (±U, 0), H = dH/dy = on y = ±L. 
Scaling these equations according to (|3.7j) with fid/dt replaced by Y/3, we arrive at 

V E = 0, ^ v4 ^ = V-(E- VH) + 2T\7 2 H -G~H, 

E = 4 [E + tr(E)l] , E = I (vu h + Vu[ + ^VHVH^j , 

where S = U /(ot 2 L) is the dimensionless applied shear strain. Comparing these equations 
with those in the previous section makes the analogy transparent. 



4. Low-dimensional viscous plate theory: buckling analysis 

Before we discuss the case of a sheared viscous plate, we begin with a brie f discussion of 
a viscous plate subject to pure compression (see Biot 1961 ; Ramberg] 1963 ) to clarify the 
role of surface tension and gravity. In this case the bounding walls move perpendicular 
rather than parallel to their length, however the same governing framework applies. 
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4.1. Pure compression 
The flat base state subject to compression is given by 

Uhb = (0, —y), H b = wb = 0, 

with compressive in-plane stresses 



-4 
-8 



where subscript b denotes base. Onset of the instability can be gleaned from the evolution 
equations for the out-of-plane displacement (|3.8b .c). which become 

for infinitesimal perturbations. We look for short-wavelength solutions and thus ignore 
lateral boundary conditions and set H (y, t) = e at cos ky. Then the growth rate a is given 
by 

*-a = 2{T-A)/e-G-/k\ 

for wavenumber k. If T > 4. then perturbations of all wavenumbers decay. Conversely, if 
T < 4, then wavenumbers k > \/G~ /2(4 — T) are unstable, with longer waves suppressed 
by buoyancy. The most unstable wavenumber is k = y / G~/(4 — T). Because onset is 
predicted to occur at infinitesimal wavelengths, the condition T = 4 may only be accurate 
for arbitrarily thin plates. However, for G~/(4 — V) not too large, the model is more 
generally valid and so the structure and wavelength of the most unstable mode should 
be captured correctly. 

4.2. Simple shear 

Shear is associated with motion parallel to the boundary and implies that we must now 
consider the evolution of infinitesimal perturbations to the flat base state 

u h b = (y,0), H b = w b = Q, (4.1a,6,c) 

having in-plane stresses 

' 2 
2 



(4.1d) 



Introducing normal mode perturbations Sf'(y)e lkx+(Tt to each variable f(x, y, t), where 
a is the growth rate, k the wavenumber and S <C 1 the amplitude, substituting into the 
governing equations (|3.8[) and linearising about the base state, we obtain the eigenvalue 
problem 

1 _ 2k 2 -^ + w' = Aik^$- +2T\ ZX- - k 2 H') - G~H' , w' = aH\ 



QyA Qy2 j dy \ dy 2 

(4.2a, &) 

subject to the boundary conditions H' = dH' /dy = on y = ±1. To find a for a given 
fc, we discretize the above equation in y using a Chebyshev pseudo-spectral method 
( Trefethen|[2000l ) and solve the resulting generalised eigenvalue problem using the eig 



routine of Matlab. 

It can be shown that (| 4 . 2 [) is self-adjoint (e.g., following an analysis similar to that 
of ISouthwell fc Skanl[l92l 7 thus a € R and all modes are stationary. However as noted 



earlier, the viscous plate model does not preserve material-frame invariance and modes 
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Figure 2. Plate model (|4,2|l . (a) Growth rates for the four most unstable modes with 
r = G~ = 0. Contours of centre-plane displacements are shown in (b), (c) for the two dominant 
modes' most unstable wavenumber. Grey/black curves indicate deflections of opposite sign. Per- 
turbations are normalised to have maximum amplitude unity and contours are equally spaced 
at intervals of 0.2. 

would also appear stationary in reference frames fixed with either wall. Again the reason 
is that in-plane advection of wrinkles contributes negligibly to the evolution of the centre- 
planc deflection compared to the out-of-plane velocity. Physically we expect that modes 
spanning the width of the plate cannot be biased by either bounding wall and thus are 
stationary in the frame of the centre-line. This is captured by reintroducing advection at 
leading order and is shown to be correct in the full Stokes calculations of iJH 

4.2.1. G~ = T = 

The only stabilising mechanisms in the viscous plate model are surface tension and 
buoyancy; without them the sheet buckles at any shear rate and all wavenumbers are 
unstable as shown in figure [2^,. Shear preferentially couples to the shortest wavelengths, 
however these waves are also most inhibited by bending resistance. Thus there is a most 
unstable mode at an intermediate wavelength, A = 27r/fc = 3.32. This mode spans the 
width of the plate and has crests and troughs aligned at roughly 45° (figure Wp)- There 
is also a cascade of subdominant modes having smaller growth rates. These differ from 
the dominant mode by having multiple crests and troughs across the width of the plate 
as shown in figure [2p. 

It is useful to compare this behaviour with the classical calculation bv lSouthwell fc Skan 
( 1924 ) for the buckling of a sheared elastic plate. The flat base state remains the same 
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Figure 3. Plate model (|4.2[) . Most unstable wavenumber (solid curves) and smallest unstable 
wavenumber (dashed curves) for varying surface tension V. Curves are shown for gravity numbers 
G~ increasing in the direction of the arrow from to 5 in intervals of unity. 

as (|4.ip (modulo the interpretation that for the elastic case, we consider displacements 
rather than velocities, and strains rather than strain rates). In contrast to the viscous 
plate, this state is stable below a non-zero threshold shear. At the onset of buckling, any 
infinitesimal out-of-plane deflection proportional to H'(y)e lkx satisfies 



with H' = dH' / dy — at y = ±1, yielding an eigenvalue problem for the critical shear 
strain. This equation is identical to that for the viscous plate on equating the reciprocal 
shear 1/S with a. In consequence, the elastic mode observed at smallest S is identical in 
structure and wavelength to the fastest growing viscous mode. 



Surface tension and buoyancy both stabilise long waves, but neither provides a short- 
wave cut-off (see figure [3]). As T and G~ increase, the long- wave cut-off shifts to increas- 
ingly short waves and in the limit r f 1, only the shortest waves remain unstable. As for 
pure compression, this is inaccurate for all but the thinnest plates because the viscous 
plate model is generally not applicable in the short-wave limit. However above onset the 
most unstable mode has moderate wavelength and thus should be faithfully reproduced. 

5. Stokes formulation and buckling analysis 

To find onset conditions for non-vanishingly thin sheets, investigate the behaviour of 
short-waves and to verify the predictions of the asymptotic viscous-plate theory of the 
previous sections, we turn to a linear stability analysis for the full Stokes equations. 

5.1. Governing equations 

The equations for conservation of mass and momentum in an incompressible fluid are 
given by 




4.2.2. G~ , T ^ 



V ■ u = 



= V • <t — pge 



(5.1a, &) 
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where the gradient operator is now three-dimensional, u = (u,v,w) is the full three- 
dimensional velocity and <r is the Cauchy stress given by 

& = -p\ + (i (Vu + Vu T ) , (5.2) 

with p pressure, superscript T denoting transpose and the gradient operator again three- 
dimensional. 

On the solid bounding walls we prescribe no-slip and no-penctration conditions 

u = (±C/,0,0) on t/ = ii, (5.3) 

while on the free surfaces at z — £ (x,y,i), we apply traction boundary conditions: on 
the upper surface, the plate is subject only to surface tension and 

6- • n + = -7 + K + n + on z = ( + 7 (5.4) 

with normal n + = (—d( + /dx,—d( + /dy,l) and curvature k + = V ■ n + (to simplify 
the presentation we do not use unit surface normals; this is permissible because we 
only consider infinitesimal deformations). On the lower surface, the plate experiences a 
pressure from the underlying fluid in addition to surface tension so that 

& ■ n = 7~K~n~ — fi~nT on z = (~, (5-5) 

with variables defined as above. In the underlying fluid, we neglect inertial and viscous 
contributions. Thus the pressure is hydrostatic and follows the relation 

-p~ = -2pgd+p~g(d+C). (5.6) 

We additionally have the kinematic conditions 

-7TF + u ~^~ + v-^r = w onz = Q , (5.7) 
dt ox dy 

where t is time. 

5.1.1. Scaled equations 

To aid comparison with results from the low-dimensional viscous plate model, we scale 
the basic variables following (|3.7p (also see Appendix [SJ, and thus set 

u = U(u,v,w/a), x = L(x,y,az), i=(a 2 L/U)t, ( ± =aL(^ ± , p=(p,U/L)p, 
(vxx,o- X y,d-yy,d- xz ,d- yz ,a zz ) = (p,U/L)(a xx i &xy i ®yy •> ™&xz •> &&yz •> ^ O ' zz 

). (5.8) 

The governing equations then become 

g du ^ dv ^ 1 dw (5 9) 

dx dy a 2 dz ' 
= V er -Ge z , (5.96) 
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with free-surface boundary conditions 
<r • n+ = -r+K+n+ 

ern = r-KTfiT - [2G- (1 + C~)G~] n - 



U> = — r h a 2 | U— h W 



on z = C + , 


(5.9cf) 


onz = f, 


(5.9e) 


on z = ( f , 


(5.9/) 







at \ dx dy 

and boundary conditions on the wall 

u= (±1,0,0) onj/ = ±l, 

where n± = (-d^/dx, -d^/dy, 1), n± = (-oPdC^/dx, -oPdC^/dy, 1), = -V 2 ^ 
and the parameters are as defined in 

5.2. Base state and perturbation equations 

We again assume a flat base state with a uni-directional, steady velocity profile and 
hydrostatic pressure: 

C 5 ± =±l, u b = (y,0,0), - Pb = a 2 G(z-l). 

To consider the evolution of infinitesimal perturbations to this base state, we assume a 
normal mode decomposition, with each dependent variable f(x,y,z,t) perturbed by an 
amount 5f'(y,z)e tkx+ ' 7t , where k is a wavenumbcr, a the growth rate and 5 < 1 the 
amplitude. Making the appropriate substitutions into the governing equations (|5.9|) and 
linearising about the base state, we obtain the eigenvalue problem for a 

iku' +^- + \^- =0, -ikp' + hu' = -^-+Lv' = -^-+Lw' = 0, (5.10o-d) 
oy or oz oy oz 

where L = — k 2 + d 2 /dy 2 + (l/a 2 )d 2 /dz 2 , subject to 





dv! , aSC^ 

— h jkuj = a — — 

az ay 


on z = 


±1, 


(5.10e) 




dv' dw' 2.7 >/+ 

«- + -»- = " * fc C 

az ay 


on z = 


±1, 


(5.10/) 


-^p'-\ 

cr 


-2^/' =r+v 2 C'+ GC' + 

a 4 oz 


on z = 


= 1, 


(5.10 ff ) 


1 . 1 dw' 
— p' + 2 — — 

or or oz 


= -r-v 2 ('- + (G~ - G)C'" 


on z = 


-1, 


(5.10/i) 




aC'* + o?ikyC,' ± =to / , 


on z = 


±1, 


(5.10i) 



and 

u = v = w = f /=t =0 on y = ±1. (5.10j-m) 

We note that as a — > this system at leading order is identical to the eigenvalue 
problem (|4.2[) for the low-dimensional viscous plate theory. In other words, the pertur- 
bation expansion and the asymptotic analysis commute. This can be shown most readily 
by starting with the system above in a stress formulation and following an equivalent 
asymptotic analysis to that given in Appendix |A] 

5.2.1. Numerical method and spectrum 

We find the eigenvalues and eigenmodes numerically, using a Chebyshev pseudospectral 
discretization in y and z and solv ing the resulting generalised linear eigenvalue problem 
using the eig routine of Matlab (jTrefethenl 120001 ) . Pressure and velocity are collocated 
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Figure 4. Stokes equations (|5.10j) . Eigenvalue spectrum for k = 8, a = 0.02, r + = T~ = 0.5 
and G = G~ — 0. Diamonds are resolved modes in the discrete spectrum, circles, o are modes 
in the continuous spectrum and crosses x are incompletely resolved modes (they still move 
perceptibly with increasing resolution). Resolution: 32 grid points in y and 16 in z. 



on the same grid points, thus to have sufficient equations we augment the boundary 
conditions on the free surfaces by the continuity equation and the bound ary conditions 
on the walls by the normal compone nt of the momentum equation ( c/. Greshol Il99ll ; 



Canuto. Hussaini fc Quarteronil 120071 ). We treat the corners as part of the walls. We 



discuss these choices further below but first describe the structure of the eigenvalue 
spectrum. 

A sample spectrum for a particular set of parameters is shown in figure |4] It consists 
of both a discrete part (indicated by dia monds) and a continuous par t (indicated by cir- 
cles), as expected for shear instabilities (jSchmid fc Henningsonll200l[ ). The discrete part 
is further broken up into purely real eigenvalues and complex-conjugate pairs. The for- 
mer correspond to stationary, sheet-spanning modes; example displacement profiles are 
shown in figures [5p, ci-iii. The latter correspond to a pair of travelling modes symmet- 
ric to one another under a 180° rotation about the vertical axis; example displacement 
profiles are shown in figures [5]C, civ, v. The continuous part of the spectrum corresponds 
to perturbations localised at a given cross-stream location y and travelling at the lo- 
cal base-state velocity. These modes are stable (although for weak stabilisation, only 
marginally so). A 'balloon' of incompletely resolved modes surrounds this continuous 
spectrum (progressively collapsing onto it with increasing resolution) and causes some 
numerical difficulties finding cut-offs where the continuous spectrum is marginally stable. 
This balloon appeared to be closest to the continuous spectrum at a given resolution for 
the augmented boundary conditions that we used, motivating our choice. 

To understand the effect of the choice of augmented boundary conditions on our results, 
we tried several different combinations of these conditions (continuity and the normal- 
component of the momentum equation) and the corner treatments (whether part of 
the wall or part of the free surface). With increasing numbers of grid points in each 
direction, the discrete spectrum visually converged to the same values for all choices. 
Similarly, using the reduced governing equations obtained by eliminating u! and p' from 
(|5.10[) also yielded the same discrete spectrum. Nevertheless, some details of the solutions 
and spectra are impacted by the choice of discretization procedure on the boundaries. 
In particular, the pressure singularity at the corner is sensitive to the treatment of the 
corner, however different combinations of the wall and free-surface boundary conditions 
only modified the solution in the grid points nearest the corners, with the remaining grid 
points' values visually appearing unaffected. The continuous spectrum is also somewhat 
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sensitive to the choice of augmented boundary conditions: specifically the size of the 
balloon varied and spurious checkerboard pressure modes appeared in the spectrum for 
continuity conditions on all boundaries. 

All solutions are presented for 28 grid-points in y and 14 points in z (unless otherwise 
stated). We calculated all eigenvalues and mode structures at this resolution and verified 
them by comparison with their counterparts at a resolution of 40 x 20. 

5.3. Dispersion relations 

To make the investigation of parameter space more manageable, we present results for 
plates stabilised only by an equal surface tension on the upper and lower surfaces, setting 
G = G~ = and F + = T~. Representative dispersion relations for the four most 
unstable modes are shown in figure [5] together with selected mid-plane displacement 
profiles. The structure of the dispersion relation is generally similar to the viscous plate 
prediction: shear couples most strongly to the shortest waves, however such modes are 
also most damped by surface tension and viscous resistance, thus the system is most 
unstable at intermediate wavenumbcr. Both long-wave and short-wave cut-offs exist. 
The longest waves are stabilised by surface tension, while the shortest are stabilised by 
an internal viscous resistance induced by shear deformations through the thickness. This 
is not accounted for in the low-dimensional plate model, which has only has extensional 
deformations through the thickness. 

For small and moderate wavenumber, the four most unstable modes are stationary 
(Imcr = in Figure [5^). However at a critical wavenumber k cr n ~ 7.5, the two most 
unstable modes have equal growth rates and for larger k they form a pair of travelling 
modes with complex conjugate growth rates. At a larger wavenumber still, the third and 
fourth modes undergo a similar bifurcation. The evolution of the mode structures with 
k reflects this changing behaviour: for small and moderate wavenumber, the dominant 
mode has a single crest or trough that spans the width of the sheet and is aligned at 
approximately 45° (figures [5]Ci-iii) . The first sub-dominant mode initially has two ex- 
trema across the sheet with a weakly deformed centre-line (figures [5j:i,ii). As k increases, 
the mid-surface deformation also increases, and eventually the mode consists of a single, 
somewhat sinuous crest or trough spanning the sheet (figure [5j:iii). At fc C rit, the first two 
modes become identical, and when k > fc cr it, they are related to one another by a 180° 
rotation about the z-axis (figures [5p, civ, v). With increasing k, they become increasingly 
concentrated on one side of the plate and travel with a fraction of the corresponding wall 
speed (figure EK). 

We see that for small and moderate wavenumbers, the viscous plate model approxi- 
mates the growth rate and mode structures of the full Stokes equations very well (bold 
grey curves in figure [5)d and final pair of rows of profiles, figures [5jC, c). However the 
transition from stationary to travelling waves at k cr n is not captured, and indeed cannot 
be captured by a regular asymptotic expansion of any order in a. 

To capture this transition requires us to reintroduce advection in the low-dimensional 
viscous plate theory. We do this by replacing (|4.2b ) with the full kinematic condition, 
which in the spectral setting reads 

w' = oil' + a 2 ikyH'. 

Predictions using this modification are included in figure [5] with short-dashed curves and 
form the central pair of rows of profiles. They provide an accurate approximation of the 
behaviour at and around the critical wavenumber as well as improving the approximation 
at smaller k. Unfortunately, this modified model also does not appear to have a short- 
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Figure 5. Stokes equations (|5.10[) . Upper panels show growth rates (a) and phase speeds (b) 
for the four most unstable modes with capillary numbers T + — T~ — 0.3, gravity numbers 
G = G~ — and aspect ratio a = 0.02. Black solid curves are the full numerical solutions, bold 
grey curves are the viscous plate prediction, black short-dashed curves are the advection-aug- 
mented plate model and black long-dashed curves are the maximum possible growth rates and 
phase speeds using the short-wave approximation. Curves are truncated where they become un- 
reliable. Lower panels show out-of-plane displacement contours at the wavenumbers indicated. 
Top pair (C,c) are full solutions, bottom (E,e) are viscous plate and middle (D,d) are the ad- 
vection-augmented plate. First of each pair (C,D,E) is the dominant mode and second (c,d,e) is 
the next most unstable. Grey/black curves indicate deflections of opposite sign. Perturbations 
are normalised to have maximum amplitude unity and contours are equally spaced in intervals 
of 0.25 (the contour is omitted for \y\ > 0.6 in C,cv). 




Figure 6. Stokes equations (|5.10|l . Growth rates for various aspect ratios a = 0.01, 0.02, 0.04, 
0.05, 0.08 (in direction of arrow) for T + = T~ = 0.6, G = G~ = 0. Black solid curves are the full 
numerical solutions, bold grey curve is the viscous plate prediction and black dashed curves are 
the advection-augmented plate model. Inset shows the critical wavenumber fe cr it and short-wave 
cut-off as functions of the aspect ratio. 



wave cut-off. However the spectrum has a neutrally stable continuous part and so it is 
difficult to state this definitively. 

Just as the low-dimensional plate theory allows us to describe moderate to long wave- 
length deformations (relative to the thickness) , one can ask if there is another asymptotic 
approximation that allows us to consider sh ort-wave behaviour (relative to the thickness) . 
The analysis of Benjamin fc Mullml (1988) provides such a route and is described in Ap- 



pendix [Bj However the predictions provide only a fair approximation to the short-wave 
behaviour (see figure [5t>) ; the discrepancy is presumably because the dominant mode 
begins to concentrate near the outer walls and so lateral boundary conditions remain 
important. 

5.4. Comparing the Stokes and asymptotic theories: dependence on aspect ratio a 

To understand how well/poorly the low-dimensional viscous plate theory does in pre- 
dicting onset of the buckling instability, we show a comparison of growth rates for the 
dominant mode at various aspect ratios in figure |6] The growth rates of the longest and 
most unstable wavelengths remain identical for aspect ratios up to 0.04 and are well de- 
scribed by the plate model (grey curve). The advection-augmented plate model (dashed 
curves) somewhat over predicts the maximum growth rates for a > 0.04, but accurately 
predicts fc cr j t and improves the estimate of the most unstable wavenumber. We can es- 
timate how fc cr jt scales with a by balancing the bending-resistance-modulated advection 
term of order k 5 a 2 in the advection-augmented model with the shear term of order k. The 
resulting prediction k 2 [it = 0(1/ a) is in good agreement with numerical solutions (inset 
of figure IBJ). Physically this results from progressively thinner plates buckling quicker, and 
thus only advection of increasingly short waves influencing mode structure. Short waves 
are cut off for fc cut = 0(l/a), when wavelengths approach the plate thickness (inset of 
figure [6|) . 
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Figure 7. Stokes equations (|5.10[) . Full numerical region of instability (shaded, bounded by 
solid black curves), most unstable wavenumber (solid black curves) and travelling waves (la- 
belled and bounded by solid black curve) as functions of F for (a) a — 0.01 and (b) 0.08 with 
G = G~ = and T + = r _ . Bold grey curves give the corresponding viscous plate predictions, 
black short-dashed curves are the advection-augmented plate model and black long-dashed are 
the short-wave approximation for the short-wave cut-off. In (a), the advection-augmented plate 
prediction for the long-wave cut-off is indistinguishable from the viscous plate. The viscous 
plate model is unable to predict the transition to travelling waves; neither plate description can 
capture the short-wave cut-off. (c) Maximum F for instability as a function of a. 

6. Stability diagrams for the onset of shear-induced buckling 

With the analysis of the different regimes at hand, we can now present the stability 
diagrams for the buckling instability of a sheared plate in terms of the dimensionless pa- 
rameters that characterise the forcing, the scaled surface tension relative to shear rate T 
and the geometrical aspect ratio a. Figure [7] shows the principal features of the instability 
for different T and a. Instability is possible for sheets of surprisingly large aspect ratio, up 
to about 0.6 (figurc[7J:). In general, the long-wave cut-off and most-unstable wavelength 
are well- approximated by the low dimensional plate model, even for reasonably thick 
plates, and the transition to travelling waves is captured by the advection-augmented 
version of the model. The short-wave cut-off is not well-described by any of the ap- 
proximations and appears to require full numerical calculation. However, portions of the 
short-wave cut-off cannot be adequately resolved even with the full numerics, because 
the most unstable mode is increasingly localised to the outer wall and thus influenced 
by the pressure singularity. This is particularly true when the continuous spectrum is 
marginally stable (r small or zero). 

7. Discussion and conclusions 

We have presented conditions for the linear, shear-induced buckling of a viscous plate 
stabilised by internal viscous resistance, surface tension and buoyancy. In the limit of a 
vanishingly thin plate, a low-dimensional asymptotic theory of the dynamics yields the 
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result that the onset of instability occurs at a scaled inverse shear rate T = 1 independent 
of buoyancy. For plates with a finite thickness stabilised only by an equal surface tension 
at the upper and lower surfaces, our numerical solutions of the eigenvalue problem based 
on the Stokes equations show that the most unstable mode has moderate wavelength: 
shear couples primarily to the shortest waves, however these are also most strongly 
suppressed by viscous resistance to bending and surface tension. This mode is stationary 
in the frame of the centre-line of the plate, spans the width of the plate, has crests aligned 
at approximately 45° and is closely related to elastic shear modes that have been well 
known for more than 85 years. 

We have compared results using the full Stokes equations and its two limiting the- 
ories, the long wav e length viscous plate model and the short-wave approximation of 
Benjamin fe Mullinl (|l988l ). The plate model predicts onset at V = 1 with infinitesimal 



wavelengths, which is only consistent and accurate for arbitrarily thin plates. However 
for plates up to a ~ 0.04, given onset parameter values, the model accurately reproduces 
both the modal structure, and the most unstable mode above onset. The short-wave 
model is somewhat inaccurate for both onset (T = 1 independent of a at an inconsistent, 
infinite wavelength) and growth rates, at least in the absence of buoyancy. 

Thus we see that the viscous plate model is quite useful, even outside its range of valid- 
ity. However, it has two short-comings. The first is that it is not material-frame invariant. 
This can be remedied by including advection at leading order, a modification which also 
extends the model's predictive power to reliably capture travelling waves that appear at 
moderately short wavelengths, even if this is not strictly correct in an asymptotic set- 
ting. The second short-coming is that amplitude-saturated modes cannot be described 
because the equations become linear when the flow is steady, and thus the amplitude is 
indeterminate. This suggests that saturated modes have larger amplitude than assumed 
in the derivat i on of t he governin g equations, and moderate or large curvature descriptions 
( HowelllflQQi 119961 : lR.ibell200lh may be necessary to describe their structure. 

Perhaps surprisingly, there are no experiments in this rectangular Couette geometry, 
even though it is close enough to many industrial flow settings associated with the float- 
glass and polymer manufacturing industries. The annular geometry, in contrast has been 
studied, although for the reasons alluded to in the introduction, that problem is funda- 
mentally different du e to th e presence of an additional length scale. Some experiments by 



Suleiman fe Munsonl (|198lh in an annular geometry do approach the rectangular limit. 



The modes observed filled the width of the annulus, were stationary with respect to the 
centre-line of the sheet, aligned at 45° to the bounding walls, and closely resembled elas- 
tic modes. Unfortunately, there are no reported data for the parameter values associated 
with the onset of the instability, so that the next step is clearly an experimental study 
of the onset of buckling in a long, rectangular Couette geometry. 



Appendix A. Long-wavelength approximation for a viscous plate 

Sever al derivations exist o f the v iscous plate model, from the analogy w ith the elastic 
plate bv lBenjamin fc Mullinl ( 1988 ) to the formal asymptotic expansion of lHowelll (|l994l 
1996h . Our attempt is to provide a more physically motivated asymptotic derivation. 



To avoid repetition, we begin with the governing Stokes equations given in ()5.1[) - (|5.7[) 
and non-dimcnsionalize them according to (|5.8[) . We justify the scalings inherent in (|5.8[) 
for the limit of the aspect ratio a —> as follows. First, in-plane stresses &h (where 
subscript h denotes in-plane components x and y) are driven by boundary motions and 
scale as fiU/L. Then for small deflections of the plate of order a, these in-plane stresses 
generate out-of-plane stresses on a cross-section a xz and a yz of order a/iU/L courtesy of 
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the in-plane components of the force balance equations and the traction conditions. These 
components in turn generate an out-of-plane stress on the surface a zz of order a 2 ^U/L 
courtesy of the out-of-plane component of the force balance equations and the traction. 
Second, there must be internal viscous resistance to deformation, specifically elongation 
on outer surfaces of wrinkles is resisted as is compression on inner surfaces, thus du/dz ~ 
dw/dx and dv/dz ~ dw/dy in the shear stresses a xz and a yz . In consequence, the out- 
of-plane velocity w appears order l/a larger than the in-plane, reflecting the fact that 
out-of-plane deformation is much easier than in-plane for a thin geometry. Thus we again 
arrive at the non-dimensional equations (|5.9[) . 

We now proceed to a regular asymptotic analysis, expanding each variable / as /o + 
oi 2 h + 0(a 4 ), substituting these into the governing equations and equating terms at each 
order in a. 

For the scalings to be adhered to in the continuity equation and the expression for a zz , 
we immediately find that the out-of-plane velocity is uniform across the plate 

dwo 
dz 

and 

= 2 ( + — 

\ dx dy 

Then the in-plane stresses become 

a h0 = 2[e + tr(e)l] , 

with e = (Vufto + ^ 7u Io)/^' * ne in-plane part of the strain rate. The expressions for 
stresses a xz and a yz imply 

u/ l0 = -zVw +Uh(x,y), 

where ii/j is the velocity on the centre-plane. The kinematic conditions on the free surfaces 
furnish 

dCo + _ dCo" 
dt dt 

Thus the plate retains uniform thickness 

C + -Co =2 

and 

w Ql (Al) 



w . 



dt 

where i/o = (Co + +Co _ )/2- 

Now evaluating forces on a cross-section, we find 

1= / o-w,d* = 4[E + tr(E)l], (A 2) 

with 

E = - (Vuft + Vul + VHVw + VwVH) , (A 3) 

the in-plane deformation rate averaged across the thickness of the plate. Integrating 
the in-plane momentum equations across the cross-section and applying the traction 
conditions, we obtain the in-plane balance of forces 



VI = 0. (A 4) 
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Similarly evaluating moments, we find 

K + 4 
M = / zo- ho dz = -- [Wwq +tr(VVw )l] + H T. 

Integrating the out-of-plane momentum equation by parts across the cross-section and 
applying boundary conditions, we obtain the out-of-plane evolution equation 

VV : M + 2rV 2 iJ - G~ Hq = 0, (A 5) 

where the term VV : M expands to -(8/3)V 4 w + V • (E • Vff )- 

The system (|A 1|) (|A 5| provides the non-dimensional viscous plate equations (|3.8[) . 
The lateral boundary conditions are obtained on setting \iho and u>o to the prescribed 
values on the walls; thus is prescribed and the normal component of Vu>o is zero. 



Appendix B. Short-wavelength approximation for a viscous slab 

In the short-wave approximation for the Stokes eigenvalue problem, we assume it is 
reasonable to ignore the lateral boundary conditions and thus we Fourier decompose in 

y as f'(y,z) = f'(z)e tmy for cross-sheet wavenumber m. Then the eigenvalue problem 
(|5.10[) reduces to 



(a/a 2 +iky) z = A 





where H' = (C' + + C' )/2 is the centre-plane deformation, h! = (C' + — C' )/2 is the 
thickness variation and 

(_ Akrn+c 2 (G~ +2K 2 T) _ 2 2G-G' + K 2 (F+ ~r~) \ 

a iK{sc-aK) aC iK(sc-aK) \ 

2 2G-G-+K 2 (r+-r~) ikva-s 2 {G-+2K i T) I 

" a iK{sc+uK) " iK{sc+aK) / 

where K = \Jk 2 + m 2 , s = sinhfqJC) and c = co sh(qJC). (These are equivalent to 



equations (27) and (28) of Benjamin k, Mullinl (|1988l ).) For our example in figure[Sl this 
translates into sinuous modes having maximum growth rate 

3 2km + c 2 K 2 T 

Re a = —a' max — — , 

m 2K(sc — aK ) 

and varicose modes, which are much less unstable (none of our modes shown in figure [5] 
are varicose), having maximum growth rate 

o -2km + s 2 K 2 T 
Re o~ = —a max ■ 



2K(sc + aK) 



The former have crests aligning with the shear, at roughly 45°, while the latter are 
anti-aligned. 

We include predictions using this approximation in figures [5] and [7] 
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